home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / answrbok / 6_11.lha / 6_11 / 6_11_div.c < prev    next >
Text File  |  1993-08-08  |  8KB  |  393 lines
  1. * Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
  2. * The C++ Answer Book */
  3. * Tony Hansen */
  4. * All rights reserved. */
  5. *
  6.    Divide u[1..m+n] by v[1..n] to form
  7.    q[0..m] and r[1..n]
  8.  
  9.    Based on:
  10.  
  11.    The Art of Computer Programming, volume 2
  12.    D. Knuth, Section 4.3.1, Algorithm D and
  13.    exercise 16.
  14. /
  15. include <arbint.h>
  16. include <debug.h>    /* DELETE */
  17.  
  18. oid dodivmod(const arbint &tmp_u, const arbint &tmp_v,
  19.    arbint "ient, arbint &remainder)
  20.  
  21.    //    Make u (the dividend) and
  22.    //    v (the divisor) positive.
  23.    int negu = tmp_u.isneg();
  24.    int negv = tmp_v.isneg();
  25.    arbint u(+tmp_u);
  26.    if (negu) u = -tmp_u;
  27.    arbint v(+tmp_v);
  28.    if (negv) v = -tmp_v;
  29. f (debug&16)    outputarb(cerr, "u=", u.p->value, u.p->length, u.p->refcnt);    // DELETE
  30. f (debug&16)    outputarb(cerr, "v=", v.p->value, v.p->length, v.p->refcnt);    // DELETE
  31.  
  32.    /*
  33. The sign of the remainder will
  34. match the sign of the dividend.
  35. The sign of the quotient will be
  36. negative if the sign of the divisor
  37. and dividend do not match, else
  38. positive.
  39.    */
  40.    int negremainder = negu;
  41.    int negquotient = (negu != negv);
  42. f (debug&16)    cerr << "neg remainder=" << negremainder << ", neg quotient=" << negquotient << "\n"; // DELETE
  43.  
  44.    // Set local variables.
  45.    ARB_type *uv = u.p->value;
  46.    ARB_type *vv = v.p->value;
  47.    int m_n = u.p->length;    // m + n
  48.    int n = v.p->length;
  49.    int m = m_n - n;
  50. f (debug&16)    cerr << "m+n=" << m_n << ", n= " << n << ", m=" << m << "\n";    // DELETE
  51.  
  52.    /*
  53. For n == 1, use simpler algorithm
  54.    */
  55.    if (n == 1)
  56. {
  57. f (debug&16)    cerr << "n==1, simpler algorithm\n";    // DELETE
  58. if (vv[0] == 0)
  59.     {
  60.     error("division by zero!");
  61.     quotient = remainder = +u;
  62.     }
  63.  
  64. else
  65.     {
  66.     /*
  67.     See Exercise 16 after Section 4.3.1
  68.     */
  69.     ARB_type *qv = new ARB_type[m_n];
  70.     ARB_Ltype prevu = 0;
  71.     ARB_type v1 = vv[0];
  72.  
  73.     for (int r = 0; r < m_n; r++)
  74.     {
  75.     ARB_Ltype t = uv[r] + prevu * ARB_base;
  76.     ARB_type tmpq = t / v1;
  77.     qv[r] = tmpq;    // % ARB_base
  78.     prevu = t - v1 * tmpq;
  79.     }
  80.  
  81.     arbint q(qv, m_n);
  82.     if (negquotient)
  83.     quotient = -q;
  84.     else
  85.     quotient = q;
  86.  
  87.     if (negremainder)
  88.     remainder = -prevu;
  89.     else
  90.     remainder = prevu;
  91.     }
  92. return;
  93. }
  94.  
  95.    /*
  96. Degenerate case of length(u) < length(v)
  97. i.e., m < 0, implying that u < v
  98.    */
  99.    else if (m_n < n)
  100. {
  101. f (debug&16)    cerr << "m_n < n, quotient <- 0\n";    // DELETE
  102. quotient = 0L;
  103. if (negremainder)
  104.     remainder = -u;
  105. else
  106.     remainder = +u;
  107. return;
  108. }
  109.  
  110.    /*
  111. Degenerate case of length(u) == length(v)
  112. i.e., m == 0, possibly implying that
  113. u < v or u == v
  114.    */
  115.    else if (m_n == n)
  116. {d
  117. bug&16)    cerr << "m_n == n\n";    // DELETE
  118. int cmp = arb_cmp(uv, vv, m_n);
  119. if (cmp < 0)    // u < v
  120.     {d
  121. bug&16)    cerr << "\tquotient == 0\n";    // DELETE
  122.     quotient = 0L;
  123.     if (negremainder)
  124.     remainder = -u;
  125.     else
  126.     remainder = +u;
  127.     return;
  128.     }
  129.  
  130. else if (cmp == 0)    // u == v
  131.     {
  132. f (debug&16)    cerr << "\tremainder == 0\n";    // DELETE
  133.     if (negquotient)
  134.     quotient = -1L;
  135.     else
  136.     quotient = 1L;
  137.     remainder = 0L;
  138.     return;
  139.     }
  140. }
  141.  
  142.    /*
  143. Now call out all of the guns from Knuth
  144.    */
  145. f (debug&16)    cerr << "all the guns\n";    // DELETE
  146.    // In rare circumstances, the first
  147.    // digit of u or v may be 0. This digit
  148.    // must now be discarded.
  149.    while ((*uv == 0) && (m_n > 1))
  150. { uv++; m_n--; }
  151.    while ((*vv == 0) && (n > 1))
  152. { vv++; n--; }
  153.    m = m_n - n;d
  154. bug&16)    cerr << "m+n=" << m_n << ", n= " << n << ", m=" << m << "\n";    // DELETE
  155. f (debug&16)    outputarb(cerr, "uv=", uv, m_n);                // DELETEd
  156. bug&16)    outputarb(cerr, "vv=", vv, n);                    // DELETE
  157.  
  158.    /*
  159. D1(a) [Normalize.]
  160.     Set d <- b/(v1+1).
  161.    */
  162.    ARB_type d = ARB_base / (vv[0] + 1);
  163. f (debug&16)    cerr << "D1: normalize\n";    // DELETEd
  164. bug&16)    cerr << "d=" << d << "\n";    // DELETEd
  165.    /*
  166. D1(b) [Normalize.]
  167.     Set (u[0]u[1]...u[m+n]) to
  168.     (u[1]u[2]...u[m+n]) * d
  169.    */
  170.    ARB_type *Ouv = uv;
  171.    ARB_type k;
  172.    uv = new ARB_type[m_n + 1];
  173.    if (d == 1)
  174. {
  175. // copy old value
  176. uv[0] = 0;
  177. memcpy((char*)&uv[1], (char*)&Ouv[0],
  178.     m_n * sizeof(ARB_type));
  179. }
  180.  
  181.    else
  182. {
  183. // multiply u by d
  184. k = 0;
  185. for (int Oi = m_n - 1, i = m_n;
  186.      i > 0;
  187.      Oi--, i--)
  188.     {
  189.     ARB_Ltype t = Ouv[Oi] * d + k;
  190.     uv[i] = t;        // % ARB_base;
  191.     k = t / ARB_base;
  192.     }
  193. uv[i] = k;
  194. }
  195. f (debug&16)    outputarb(cerr, "uv=", uv, m_n+1);    // DELETE
  196.  
  197.    /*
  198. D1(c) [Normalize.]
  199.     Set (v[1]v[2]...v[n]) to
  200.     (v[1]v[2]...v[n]) * d
  201.    */
  202.    ARB_type *Ovv = vv;
  203.    vv = new ARB_type[n];
  204.    if (d == 1)
  205. {d    // copy old value
  206. memcpy((char*)&vv[0], (char*)&Ovv[0],
  207.     n * sizeof(ARB_type));
  208. }
  209.  
  210.    else
  211. {
  212. // multiply v by d
  213. k = 0;
  214. for (int Oi = n - 1, i = n - 1;
  215.      i >= 0;
  216.      Oi--, i--)
  217.     {
  218.     ARB_Ltype t = Ovv[Oi] * d + k;
  219.     vv[i] = t;        // % ARB_base;
  220.     k = t / ARB_base;
  221.     }
  222. }d
  223. bug&16)    outputarb(cerr, "vv=", vv, n);    // DELETE
  224.  
  225.    /*
  226. D2 [Initialize j]
  227.     Set j <- 0
  228.  
  229. D7 [Loop on j]
  230.     Set j <- j+1
  231.     Loop if j <= m
  232.    */
  233.    ARB_type v1 = vv[0];
  234.    ARB_type v2 = vv[1];
  235.    ARB_type *qv = new ARB_type[m + 1];
  236.    ARB_type *nv = new ARB_type[n + 1];
  237.  
  238.    for (int j = 0; j <= m; j++)
  239. {
  240. /*
  241.     D3 [Calculate q^]
  242.     If u[j] == v[1]
  243.         Set q^ <- base - 1
  244.     Else
  245.         Set q^ <-
  246.         (u[j] * base + u[j+1]) / v[1]
  247.     While v[2] * q^ >
  248.        (u[j] * base + u[j+1] - q^ * v[1]) *
  249.         base + u[j+2]
  250.         Set q^ <- q^ - 1
  251. */
  252. ARB_Ltype q_hat;
  253.  
  254. if (uv[j] == v1)
  255.     q_hat = ARB_base - 1;
  256.  
  257. else
  258.     q_hat = (uv[j] * ARB_base + uv[j+1]) / v1;d
  259. bug&16)    cerr << "D3: calculate q^, j=" << j << "\n";    // DELETEdif (debug&16)    cerr << "\tq^=" << q_hat << "\n";        // DELETE
  260.  
  261. ARB_Ltype u_j = uv[j];
  262. ARB_Ltype u_j1 = uv[j+1];
  263. ARB_Ltype u_j2 = uv[j+2];
  264.  
  265. for ( ; ; q_hat--)
  266.     {d        /*
  267.     if ((v2 * q_hat) <=
  268.         ((u_j * ARB_base + u_j1 -
  269.          v1 * q_hat) * ARB_base + u_j2))
  270.         q^--;
  271.     */
  272.     ARB_Ltype u_j_q_hat = u_j * ARB_base + u_j1;
  273.     u_j_q_hat -= v1 * q_hat;
  274.     if (u_j_q_hat / ARB_base != 0)
  275.     break;
  276.  
  277.     u_j_q_hat *= ARB_base;
  278.     u_j_q_hat += u_j2;
  279.     if ((v2 * q_hat) <= u_j_q_hat)
  280.     break;
  281.     }
  282. f (debug&16)    cerr << "q^=" << q_hat << "\n";            // DELETE
  283.  
  284. /*
  285.     D4 [Multiply and subtract.]
  286.     Replace u[j..j+n] by
  287.         u[j..j+n] -
  288.          q^ * v[1..n]
  289. */d
  290. bug&16)    cerr << "D4: multiply and subtract\n";    // DELETEd    // set nv <- q^ * (v[1..n])
  291. k = 0;
  292. for (int dl = n, vl = n-1; vl >= 0; vl--, dl--)
  293.     {
  294.     ARB_Ltype t = vv[vl] * q_hat + k;
  295.     nv[dl] = t;        // % ARB_base;
  296.     k = t / ARB_base;
  297.     }
  298. nv[0] = k;
  299.  
  300. // subtract nv[0..n] from u[j..j+n]
  301. int borrow = 0;
  302. int ul = j + n;
  303. for (dl = n; dl >= 0; dl--, ul--)
  304.     {
  305.     ARB_Ltype t = uv[ul] - nv[dl] - borrow;
  306.     uv[ul] = t;        // % ARB_base
  307.     borrow = (t / ARB_base) ? 1 : 0;
  308.     }
  309.  
  310. /*
  311.     D5 [Test remainder]
  312.     Set q[j] <- q^
  313.     if u[j] < 0
  314.  
  315.     D6 [Add back]
  316.         add 0v[1..n] back to u[j..j+n]
  317.         q[j]--
  318. */
  319. qv[j] = q_hat;d
  320. bug&16)    cerr << "D5: test remainder\n";    // DELETEdif (debug&16)    cerr << "\tqv[" << j << "]=" << q_hat << "\n";    // DELETEdif (debug&16)    cerr << "\tborrow=" << borrow << "\n";    // DELETE
  321. if (borrow != 0)
  322.     {
  323.     for (k = 0, ul = j + n, vl = n - 1;
  324.      vl >= 0;
  325.      vl--, ul--)
  326.     {d        ARB_Ltype t = uv[ul] + vv[vl] + k;
  327.     uv[ul] = t;    // % ARB_base
  328.     k = t / ARB_base;
  329.     }
  330.  
  331.     uv[j] += k;
  332.     qv[j]--;d
  333. bug&16)    cerr << "\tqv[" << j << "]=" << q_hat << "\n";    // DELETE
  334.     }
  335. }
  336.  
  337.    /*
  338. D8 [Unnormalize]
  339.     q[0..m] is quotient
  340.     u[m+1..m+n] / d is remainder
  341.    */
  342. f (debug&16)    outputarb(cerr, "q=", qv, m+1);    // DELETE
  343.    arbint q(qv, m+1);
  344.    if (negquotient)
  345. quotient = -q;
  346.    else
  347. quotient = q;
  348. f (debug&16)    outputarb(cerr, "quotient=", quotient.p->value, quotient.p->length);    // DELETE
  349.  
  350.    // divide u[m+1..m+n] by d
  351.    ARB_type *rem = new ARB_type[n];
  352.  
  353.    if (d == 1)        // nothing special to do
  354. memcpy((char*)rem, (char*)&uv[m+1],
  355.     n * sizeof(ARB_type));
  356.  
  357.    elsRB_type))
  358. // do division by single digit
  359. for (int rl = 0, ul = m + 1;
  360.      ul <= m_n;
  361.      ul++, rl++)
  362.     {
  363.     ARB_Ltype t = uv[ul] + prevu * ARB_base;
  364.     ARB_type tmpq = t / d;
  365.     rem[rl] = tmpq;    // % ARB_base
  366.     prevu = t - d * tmpq;
  367.     }
  368. }
  369.  
  370. f (debug&16)    outputarb(cerr, "rem=", rem, n);    // DELETE
  371.    arbint r(rem, n);
  372.    if (negremainder)
  373. remainder = -r;
  374.    else
  375. remainder = r;d
  376. bug&16)    outputarb(cerr, "remainder=", remainder.p->value, remainder.p->length);    // DELETE
  377.    delete uv; delete vv; delete nv;
  378.  
  379.  
  380. rbint operator%(const arbint &u, const arbint &v)
  381.  
  382.    arbint quot, rem;
  383.    dodivmod(u, v, quot, rem);
  384.    return rem;
  385.  
  386.  
  387. rbint operator/(const arbint &u, const arbint &v)
  388.  
  389.    arbint quot, rem;
  390.    dodivmod(u, v, quot, rem);
  391.    return quot;
  392.  
  393.